
module vert_turb_driver_mod

!-----------------------------------------------------------------------
!
!       driver for compuing vertical diffusion coefficients
!
!         choose either:
!              1) mellor-yamada 2.5 (with tke)
!              2) non-local K scheme
!              3) entrainment and diagnostic turbulence (edt) from
!                 Bretherton and Grenier
!
!-----------------------------------------------------------------------
!---------------- modules ---------------------


use      my25_turb_mod, only: my25_turb_init, my25_turb_end,  &
                              my25_turb, tke_surf, get_tke,   &
                              my25_turb_restart

use       tke_turb_mod, only: tke_turb_init, tke_turb_end, tke_turb

use    diffusivity_mod, only: diffusivity, molecular_diff

use            edt_mod, only: edt_init, edt, edt_end

use physics_radiation_exch_mod, only: exchange_control_type
use  physics_types_mod, only: physics_control_type

use   shallow_conv_mod, only: shallow_conv_init, shallow_conv

use stable_bl_turb_mod, only: stable_bl_turb_init, stable_bl_turb

use        entrain_mod, only: entrain_init, entrain, entrain_end

use   diag_manager_mod, only: register_diag_field, send_data

use   time_manager_mod, only: time_type, get_time, operator(-)

use      constants_mod, only: rdgas, rvgas, kappa, grav
 
use            mpp_mod, only: input_nml_file
use            fms_mod, only: mpp_pe, mpp_root_pe, stdlog, &
                              error_mesg, open_namelist_file, file_exist, &
                              check_nml_error, close_file, FATAL, &
                              write_version_number, & 
                              stdout, mpp_chksum
 

use  field_manager_mod, only: MODEL_ATMOS

use tracer_manager_mod, only: get_tracer_index, &
                              get_number_tracers

use moist_proc_utils_mod,  only: rh_calc

use atmos_cmip_diag_mod, only: register_cmip_diag_field_2d, &
                               register_cmip_diag_field_3d, &
                               send_cmip_data_3d, &
                               cmip_diag_id_type, &
                               query_cmip_diag_id

implicit none
private

!---------------- interfaces ---------------------

public   vert_turb_driver_init, vert_turb_driver_end, vert_turb_driver
public   vert_turb_driver_restart


!-----------------------------------------------------------------------
!--------------------- version number ----------------------------------

character(len=128) :: version = '$Id$'
character(len=128) :: tagname = '$Name$'
logical            :: module_is_initialized = .false.

!-----------------------------------------------------------------------
 real, parameter :: p00    = 1000.0E2
 real, parameter :: p00inv = 1./p00
 real, parameter :: d622   = rdgas/rvgas
 real, parameter :: d378   = 1.-d622
 real, parameter :: d608   = d378/d622

!---------------- private data -------------------

 real :: gust_zi = 1000.   ! constant for computed gustiness (meters)

 integer :: nql, nqi, nqa    !  tracer indices for stratiform clouds
 integer :: ntke             !  tracer index for TKE
 integer :: ntp              !  number of prognostic tracers

 integer         :: outunit

 logical :: doing_prog_clouds
 logical :: use_tau

!-----------------------------------------------------------------------
!-------------------- namelist -----------------------------------------

 logical :: do_shallow_conv  = .false.
 logical :: do_mellor_yamada = .true.
 logical :: do_tke_turb      = .false.
 logical :: do_diffusivity         = .false.
 logical :: do_molecular_diffusion = .false.
 logical :: do_edt                 = .false.
 logical :: do_stable_bl     = .false.
 logical :: do_entrain    = .false.
 logical :: do_simple = .false. 

 character(len=24) :: gust_scheme  = 'constant' ! valid schemes are:
                                                !   => 'constant'
                                                !   => 'beljaars'
 real              :: constant_gust = 1.0
 real              :: gust_factor   = 1.0

 integer           :: alternate_zpbl = 0        ! alternate algorith to compute PBL height

!-->h1g, 2012-07-16
 integer :: do_clubb
 integer :: nwp2
 real    :: wp2_min = 4.e-4
 real    :: diff_min = 1.e-3    ! minimum value of a diffusion 
                                ! coefficient beneath which the
                                ! coefficient is reset to zero
 
 integer :: id_clubb_on
 integer :: id_stable_on                     ! ( diff_m_stab > diff_m  or diff_t_stab > diff_t)
 
 integer :: id_stable_effective              !     ( diff_m_stab > diff_m  or diff_t_stab > diff_t) 
                                             ! and ( diff_m_stab > diff_min  or diff_t_stab > diff_min)
!<--h1g, 2012-07-16 
 
 namelist /vert_turb_driver_nml/ do_shallow_conv, do_mellor_yamada, &
                                 do_tke_turb, &
                                 gust_scheme, constant_gust,          &
                                 do_molecular_diffusion, do_stable_bl, &
                                 do_diffusivity, do_edt, do_entrain, &
                                 gust_factor, do_simple, wp2_min, &
                                 alternate_zpbl   ! cjg: PBL depth mods

!-------------------- diagnostics fields -------------------------------

integer :: id_tke,    id_lscale, id_lscale_0, id_z_pbl, id_gust,  &
           id_diff_t, id_diff_m, id_diff_sc, id_z_full, id_z_half,&
           id_uwnd,   id_vwnd,   id_diff_t_stab, id_diff_m_stab,  &
           id_diff_t_entr, id_diff_m_entr,                        &
           id_z_Ri_025, id_tref, id_qref, id_rh_Ri_025  ! cjg: PBL depth mods, h1g, add RH diagnostics at Ri_025, 2015-04-02
integer :: id_bldep ! cmip6 boundary layer depth
type(cmip_diag_id_type) :: ID_evu, ID_edt

real :: missing_value = -999.

character(len=9) :: mod_name = 'vert_turb'

!-----------------------------------------------------------------------

contains

!#######################################################################

subroutine vert_turb_driver (is, js, Time, Time_next, dt, tdtlw,       &
                             frac_land,                                &
                             p_half, p_full, z_half, z_full,           &
                             t_ref, q_ref,                             &  ! cjg: PBL depth mods
                             u_star, b_star, q_star, rough,            &
                             lat, convect,                             &
                             u, v, t, q, r, um, vm, tm, qm, rm, rdiag, &
                             udt, vdt, tdt, qdt, rdt, diff_t, diff_m,  &
                             gust, z_pbl, mask, kbot, tke_avg          )  ! h1g: output averaged TKE within PBL  

!-----------------------------------------------------------------------
integer,         intent(in)         :: is, js
type(time_type), intent(in)         :: Time, Time_next
   real,         intent(in)         :: dt
   real, intent(in), dimension(:,:) :: t_ref, q_ref,  &  ! cjg: PBL depth mods
                                       frac_land, u_star, b_star, q_star, rough, lat
logical, intent(in), dimension(:,:) :: convect       
   real, intent(in), dimension(:,:,:) :: tdtlw, p_half, p_full, &
                                         z_half, z_full, &
                                         u, v, t, q, um, vm, tm, qm, &
                                         udt, vdt, tdt, qdt
   real, intent(in),    dimension(:,:,:,:) :: r, rm, rdt
   real, intent(inout), dimension(:,:,:,ntp+1:) :: rdiag
   real, intent(out),   dimension(:,:,:) :: diff_t, diff_m
   real, intent(out),   dimension(:,:)   :: gust, z_pbl 
   real, intent(in),optional, dimension(:,:,:) :: mask
integer, intent(in),optional, dimension(:,:) :: kbot

!---> h1g, 2015-08-11
  real, intent(out), optional, dimension(:,:) :: tke_avg  !averaged TKE within PBL
!<--- h1g, 2015-08-11

!-----------------------------------------------------------------------
real   , dimension(size(t,1),size(t,2),size(t,3))   :: ape, thv
logical, dimension(size(t,1),size(t,2),size(t,3)+1) :: lmask
real   , dimension(size(t,1),size(t,2),size(t,3)+1) :: el, diag3
real   , dimension(size(t,1),size(t,2),size(t,3)+1) :: tke
real   , dimension(size(t,1),size(t,2))             :: stbltop
real   , dimension(size(t,1),size(t,2))             :: z_Ri_025    ! cjg: PBL depth mods

real   , dimension(size(t,1),size(t,2))             :: RH_Ri_025   ! h1g: relative humidity at Ri_025, 2015-04-02

real   , dimension(size(t,1),size(t,2),size(t,3))   :: RH_3D_tmp   ! h1g: 3D relative humidity, 2015-04-02

real   , dimension(size(t,1),size(t,2))             :: el0, vspblcap
real   , dimension(size(diff_t,1),size(diff_t,2), &
                                  size(diff_t,3))   :: diff_sc,     &
                                                       diff_t_stab, &
                                                       diff_m_stab, &
       diff_t_entr, &
       diff_m_entr, &
       use_entr
real   , dimension(size(t,1),size(t,2),size(t,3))   :: tt, qq, uu, vv
real   , dimension(size(t,1),size(t,2),size(t,3))   :: qlin, qiin, qain
real    :: dt_tke
integer :: ie, je, nlev, sec, day, nt
integer :: i,j,kk
logical :: used
!-->h1g, 2012-08-07
real   , dimension(size(diff_t,1),size(diff_t,2), &
                                  size(diff_t,3))   :: clubb_on,   &
                                                       stable_on,  &
                                                       stable_effective
!<--h1g, 2012-08-07 

!-----------------------------------------------------------------------
!----------------------- vertical turbulence ---------------------------
!-----------------------------------------------------------------------

      if (.not. module_is_initialized)  call error_mesg  &
                     ('vert_turb_driver in vert_turb_driver_mod',  &
                      'initialization has not been called', FATAL)

     nlev = size(p_full,3)
     ie = is + size(p_full,1) - 1
     je = js + size(p_full,2) - 1

     if ( present(tke_avg) ) tke_avg = 0.0   ! h1g, 2015-08-11
!-----------------------------------------------------------------------
!---- set up state variable used by this module ----

      if (use_tau) then
      !-- variables at time tau
          uu = u
          vv = v
          tt = t
          qq = q
      else
      !-- variables at time tau+1
          uu = um + dt*udt
          vv = vm + dt*vdt
          tt = tm + dt*tdt
          qq = qm + dt*qdt
      endif

      !------ setup cloud variables: ql & qi & qa -----
      if (doing_prog_clouds) then
           nt=size(r,4)
           if (nt == 0 .or. nt < max(nql,nqi,nqa))                    &
        call error_mesg ('vert_turb_driver',                  & 
                     'number of tracers less than nql or nqi or nqa', &
      FATAL) 
           if (use_tau) then
                qlin (:,:,:)=r(:,:,:,nql)
                qiin (:,:,:)=r(:,:,:,nqi)
                qain (:,:,:)=r(:,:,:,nqa)
           else
                qlin (:,:,:)=rm(:,:,:,nql)+rdt(:,:,:,nql)*dt
                qiin (:,:,:)=rm(:,:,:,nqi)+rdt(:,:,:,nqi)*dt
                qain (:,:,:)=rm(:,:,:,nqa)+rdt(:,:,:,nqa)*dt
           endif
      else
           qlin = 0.0
           qiin = 0.0
           qain = 0.0
      end if

!--------------------------------------------------------------------

!--------------------------------------------------------------------
! initialize output

   diff_t = 0.0
   diff_m = 0.0
   el     = 0.0
   el0    = 0.0
   z_pbl = -999.0
   
!-------------------------------------------------------------------
! initiallize variables   
   vspblcap = 0.0   
   
!-----------------------------------------------------------------------
if (do_mellor_yamada) then

!    ----- virtual temp ----------
     ape(:,:,:)=(p_full(:,:,:)*p00inv)**(-kappa)
     if(do_simple) then 
       thv(:,:,:)=tt(:,:,:)*ape(:,:,:)
     else
       thv(:,:,:)=tt(:,:,:)*(qq(:,:,:)*d608+1.0)*ape(:,:,:)
     endif  
     if (present(mask)) where (mask < 0.5) thv = 200.

 endif

!---------------------------
 if (do_mellor_yamada) then
!---------------------------

!    ----- time step for prognostic tke calculation -----
     call get_time (Time_next-Time, sec, day)
     dt_tke = real(sec+day*86400)

!    --------------------- update tke-----------------------------------
!    ---- compute surface tke --------
!    ---- compute tke, master length scale (el0),  -------------
!    ---- length scale (el), and vert mix coeffs (diff_t,diff_m) ----

     call tke_surf  (is, js, u_star, kbot=kbot)



     if ( id_z_pbl > 0 .or. id_bldep > 0 ) then
     !------ compute pbl depth from k_profile if diagnostic needed -----
     call my25_turb (is, js, dt_tke, frac_land, p_half, p_full, thv, uu, vv, &
                     z_half, z_full, rough,   &
                     el0, el, diff_m, diff_t, &
                     mask=mask, kbot=kbot, &
                     ustar=u_star,bstar=b_star,h=z_pbl)
     else
     call my25_turb (is, js, dt_tke, frac_land, p_half, p_full, thv, uu, vv, &
                     z_half, z_full, rough,   &
                     el0, el, diff_m, diff_t, &
                     mask=mask, kbot=kbot)
     end if

!---------------------------
 else if (do_tke_turb) then
!---------------------------

!-->cjg debug
!100 format("BEFORE TURB:",A32," = ",Z20)
!  outunit = stdout()
! write(outunit,100) 't                ', mpp_chksum(t)
! write(outunit,100) 'q                ', mpp_chksum(q)
! write(outunit,100) 'z_full           ', mpp_chksum(z_full)
! write(outunit,100) 'z_half           ', mpp_chksum(z_half)
! write(outunit,100) 'qa               ', mpp_chksum(rdiag(:,:,:,nqa))
! write(outunit,100) 'tke              ', mpp_chksum(rdiag(:,:,:,ntke))
! write(outunit,100) 'el0              ', mpp_chksum(el0)
! write(outunit,100) 'el               ', mpp_chksum(el)
! write(outunit,100) 'diff_m           ', mpp_chksum(diff_m)
! write(outunit,100) 'diff_t           ', mpp_chksum(diff_t)
! write(outunit,100) 'z_pbl            ', mpp_chksum(z_pbl)
!<--cjg debug

!    ----- time step for prognostic tke calculation -----
     call get_time (Time_next-Time, sec, day)
     dt_tke = real(sec+day*86400)

!    --------------------- update tke-----------------------------------
!    ---- compute tke, master length scale (el0),  -------------
!    ---- length scale (el), and vert mix coeffs (diff_t,diff_m) ----

     if( present(tke_avg) ) then
      call tke_turb (is, ie, js, je, Time_next, dt_tke, frac_land,      &
                     p_half, p_full, z_half, z_full,                    &
                     tt, qq, qain, qlin, qiin, uu, vv,                  &
                     rough, u_star, b_star,                             &
                     rdiag(:,:,:,ntke),                                 &
                     el0, el, diff_m, diff_t, z_pbl, tke_avg=tke_avg)
     else
      call tke_turb (is, ie, js, je, Time_next, dt_tke, frac_land,      &
                     p_half, p_full, z_half, z_full,                    &
                     tt, qq, qain, qlin, qiin, uu, vv,                  &
                     rough, u_star, b_star,                             &
                     rdiag(:,:,:,ntke),                                 &
                     el0, el, diff_m, diff_t, z_pbl)
     endif   ! h1g, 2015-08-11
!-->cjg debug
!101 format("AFTER TURB: ",A32," = ",Z20)
!  outunit = stdout()
! write(outunit,101) 't                ', mpp_chksum(t)
! write(outunit,101) 'q                ', mpp_chksum(q)
! write(outunit,101) 'z_full           ', mpp_chksum(z_full)
! write(outunit,101) 'z_half           ', mpp_chksum(z_half)
! write(outunit,101) 'qa               ', mpp_chksum(rdiag(:,:,:,nqa))
! write(outunit,101) 'tke              ', mpp_chksum(rdiag(:,:,:,ntke))
! write(outunit,101) 'el0              ', mpp_chksum(el0)
! write(outunit,101) 'el               ', mpp_chksum(el)
! write(outunit,101) 'diff_m           ', mpp_chksum(diff_m)
! write(outunit,101) 'diff_t           ', mpp_chksum(diff_t)
! write(outunit,101) 'z_pbl            ', mpp_chksum(z_pbl)
!<--cjg debug

!---------------------------
 else if (do_diffusivity) then
!--------------------------------------------------------------------
!----------- compute molecular diffusion, if desired  ---------------

    if (do_molecular_diffusion) then
      call molecular_diff (tt, p_half, diff_m, diff_t)
    else
      diff_m = 0.0
      diff_t = 0.0
    endif

!---------------------------
!------------------- non-local K scheme --------------


    call diffusivity ( tt, qq, uu, vv, p_full, p_half, z_full, z_half,   &
                       u_star, b_star, z_pbl, diff_m, diff_t, &
                       kbot = kbot)

!---------------------------
else if (do_edt) then
!----------------------------

!    ----- time step for prognostic tke calculation -----
      call get_time (Time_next-Time, sec, day)
      dt_tke = real(sec+day*86400)
 

      tke = 0.0

    call edt(is,ie,js,je,dt_tke,Time_next,tdtlw, u_star,b_star,q_star, &
             tt,qq,  &
             qlin,qiin,qain,uu,vv,z_full,p_full,z_half,p_half,stbltop, &
             diff_m,diff_t,z_pbl,kbot=kbot,tke=tke)


 endif
 


 
!------------------------------------------------------------------
! --- boundary layer entrainment parameterization

   if( do_entrain ) then

       call entrain(is,ie,js,je,Time_next,tdtlw, convect,u_star,b_star,&
                    tt,qq, &
            qlin,qiin,qain,uu,vv,z_full,p_full,z_half,p_half,diff_m,   &
    diff_t,diff_m_entr,diff_t_entr,use_entr,z_pbl,vspblcap,    &
    kbot=kbot)
   
   endif

!-----------------------------------------------------------------------
! --- stable boundary layer parameterization

! --> h1g, 2012-08-08, pre-set stable_on is 0
   stable_on = 0.0
   stable_effective = 0.0
! <-- h1g, 2012-08-08

   if( do_stable_bl ) then

        if (do_entrain) then

CALL STABLE_BL_TURB( is, js, Time_next, tt, qq, qlin, qiin, uu,&
                     vv, z_half, z_full, u_star, b_star, lat,  &
     diff_m_stab, diff_t_stab,                 &
     vspblcap = vspblcap, kbot=kbot)
     
            diff_m = use_entr*diff_m_entr + (1-use_entr)*diff_m_stab
            diff_t = use_entr*diff_t_entr + (1-use_entr)*diff_t_stab
    
            !for diagnostic purposes only, save the stable_bl_turb
            !coefficient only where it was used
    
            diff_m_stab = (1-use_entr)*diff_m_stab
            diff_t_stab = (1-use_entr)*diff_t_stab    
         
else

CALL STABLE_BL_TURB( is, js, Time_next, tt, qq, qlin, qiin, uu,&
                     vv, z_half, z_full, u_star, b_star, lat,  &
     diff_m_stab, diff_t_stab,kbot=kbot)

! --->h1g, 2012-07-16
     if(  do_clubb > 0 ) then
        clubb_on = 1.0
        where ( rdiag(:,:,:, nwp2) <= wp2_min )
            where( diff_m_stab > diff_m .or. diff_t_stab > diff_t )
               stable_on = 1.0
               where( diff_m_stab >= diff_min .or. diff_t_stab >= diff_min )
                   stable_effective = 1.0
               endwhere
            endwhere
            diff_m = diff_m +  MAX( diff_m_stab - diff_m, 0.0 )
            diff_t = diff_t +  MAX( diff_t_stab - diff_t, 0.0 )
            clubb_on = 0.0
        endwhere
     else
        diff_m = diff_m +  MAX( diff_m_stab - diff_m, 0.0 )
        diff_t = diff_t +  MAX( diff_t_stab - diff_t, 0.0 )
     endif
! <---h1g, 2012-07-16

end if
        
    endif
   
!-----------------------------------------------------------------------
!------------------ shallow convection ???? ----------------------------

   if (do_shallow_conv) then
        call shallow_conv (tt, qq, p_full, p_half, diff_sc, kbot)
        diff_t = diff_t + diff_sc
   endif

!-----------------------------------------------------------------------
!------------- define gustiness ------------

     if ( trim(gust_scheme) == 'constant' ) then
          gust = constant_gust
     else if ( trim(gust_scheme) == 'beljaars' ) then
!    --- from Beljaars (1994) and Beljaars and Viterbo (1999) ---
          where (b_star > 0.)
             gust = gust_factor * (u_star*b_star*gust_zi)**(1./3.)
          elsewhere
             gust = 0.
          endwhere
     endif

!-----------------------------------------------------------------------
!------------------------ diagnostics section --------------------------

if (do_mellor_yamada .or. do_tke_turb) then

!     --- set up local mask for fields with surface data ---
      if ( present(mask) ) then
         lmask(:,:,1)        = .true.
         lmask(:,:,2:nlev+1) = mask(:,:,1:nlev) > 0.5
      else
         lmask = .true.
      endif

!------- tke --------------------------------
      if ( id_tke > 0 ) then
         call get_tke(is,ie,js,je,tke)
         used = send_data ( id_tke, tke, Time_next, is, js, 1, &
                            mask=lmask )
      endif

!------- length scale (at half levels) ------
      if ( id_lscale > 0 ) then
         used = send_data ( id_lscale, el, Time_next, is, js, 1,  &
                            mask=lmask )
      endif

!------- master length scale -------
      if ( id_lscale_0 > 0 ) then
         used = send_data ( id_lscale_0, el0, Time_next, is, js )
      endif

end if

if (do_edt) then 
    
!     --- set up local mask for fields with surface data ---
    if ( present(mask) ) then
          lmask(:,:,1)        = .true.
          lmask(:,:,2:nlev+1) = mask(:,:,1:nlev) > 0.5
     else   
        lmask = .true.
       endif

!------- tke --------------------------------
      if ( id_tke > 0 ) then
        used = send_data ( id_tke, tke, Time_next, is, js, 1,     &
                          mask=lmask )
      endif
 
end if

!-->cjg: addition for new PBL depth diagnostic

!------- z_Ri_025: bulk Richardson derived height  -------
      if ( id_tref > 0 ) then
         used = send_data ( id_tref, t_ref, Time_next, is, js )
      endif
      if ( id_qref > 0 ) then
         used = send_data ( id_qref, q_ref, Time_next, is, js )
      endif
      if ( id_z_Ri_025 > 0 .or. alternate_zpbl == 1) then
         call bulk_Ri_height_b(tt,qq,uu,vv,t_ref,q_ref,p_full,z_full,p_half,z_half,z_Ri_025)
         used = send_data ( id_z_Ri_025, z_Ri_025, Time_next, is, js )
         if ( alternate_zpbl == 1 ) z_pbl = z_Ri_025

!--->h1g: calculate relative humidity, 2015-04-02
         if ( id_rh_Ri_025 > 0 ) then
           call rh_calc ( p_full,  tt, qq, RH_3D_tmp, .false., do_cmip=.true.)
        
           rh_Ri_025(:,:) = missing_value

           do i = 1, size(z_full,1) 
             do j = 1, size(z_full,2)
!      Vertical upward loop
               do kk = nlev, 1, -1 
                 if ( z_full(i,j,kk) >= z_Ri_025(i,j) + z_half(i,j,nlev+1) ) then
                   rh_Ri_025(i,j) = RH_3D_tmp(i,j,kk) * 100.
                   exit
                 endif
               enddo
             enddo
           enddo
           used = send_data ( id_rh_Ri_025, rh_Ri_025, Time_next, is, js )
         endif  ! if  id_rh_Ri_025 > 0
!<---h1g,  2015-04-02
      endif
!<--cjg

!------- boundary layer depth -------
      if ( id_z_pbl > 0 ) then
         used = send_data ( id_z_pbl, z_pbl, Time_next, is, js )
      endif
!------- cmip name ------
      if ( id_bldep > 0 ) then
         used = send_data ( id_bldep, z_pbl, Time_next, is, js )
      endif

!------- gustiness -------
      if ( id_gust > 0 ) then
         used = send_data ( id_gust, gust, Time_next, is, js )
      endif


!------- output diffusion coefficients ---------

  if ( id_diff_t > 0 .or. id_diff_m > 0 .or. id_diff_sc > 0 .or. &
       id_diff_t_stab > 0 .or. id_diff_m_stab > 0 .or.           &
       id_diff_t_entr > 0 .or. id_diff_m_entr > 0 .or.           &
       query_cmip_diag_id(ID_edt) .or. query_cmip_diag_id(ID_evu) ) then
!       --- set up local mask for fields without surface data ---
        if (present(mask)) then
            lmask(:,:,1:nlev) = mask(:,:,1:nlev) > 0.5
            lmask(:,:,nlev+1) = .false.
        else
            lmask(:,:,1:nlev) = .true.
            lmask(:,:,nlev+1) = .false.
        endif
!       -- dummy data at surface --
        diag3(:,:,nlev+1)=0.0
  endif

!------- diffusion coefficient for heat/moisture -------
   if ( id_diff_t > 0 .or. query_cmip_diag_id(ID_edt) ) then
      diag3(:,:,1:nlev) = diff_t(:,:,1:nlev)
      if (id_diff_t > 0) used = send_data ( id_diff_t, diag3, Time_next, is, js, 1, mask=lmask )
      if (query_cmip_diag_id(ID_edt)) used = send_cmip_data_3d (ID_edt, diag3, Time_next, is, js, 1, mask=lmask)
   endif

!------- diffusion coefficient for momentum -------
   if ( id_diff_m > 0 .or. query_cmip_diag_id(ID_evu)) then
      diag3(:,:,1:nlev) = diff_m(:,:,1:nlev)
      if (id_diff_m > 0) used = send_data ( id_diff_m, diag3, Time_next, is, js, 1, mask=lmask )
      if (query_cmip_diag_id(ID_evu)) used = send_cmip_data_3d (ID_evu, diag3, Time_next, is, js, 1, mask=lmask)
   endif

!------- diffusion coefficient for shallow conv -------
 if (do_shallow_conv) then
   if ( id_diff_sc > 0 ) then
      diag3(:,:,1:nlev) = diff_sc(:,:,1:nlev)
      used = send_data ( id_diff_sc, diag3, Time_next, is, js, 1, mask=lmask)
   endif
 endif

!------- diffusion coefficients for stable boudary layer -------
   if (do_stable_bl) then
!------- for heat/moisture -------
    if ( id_diff_t_stab > 0 ) then
       diag3(:,:,1:nlev) = diff_t_stab(:,:,1:nlev)
      used = send_data ( id_diff_t_stab, diag3, Time_next, is, js, 1, mask=lmask )
  endif
!------- for momentum -------
    if ( id_diff_m_stab > 0 ) then
       diag3(:,:,1:nlev) = diff_m_stab(:,:,1:nlev)
     used = send_data ( id_diff_m_stab, diag3, Time_next, is, js, 1, mask=lmask )
    endif

     if ( id_stable_on > 0 ) then
        used = send_data ( id_stable_on, stable_on, Time_next, is, js, 1, mask=lmask )
     endif

     if ( id_stable_effective > 0 ) then
        used = send_data ( id_stable_effective, stable_effective, Time_next, is, js, 1, mask=lmask )
     endif

 endif

!------- diffusion coefficients for entrainment module -------
 if (do_entrain) then
      if ( id_diff_t_entr > 0 ) then
       diag3(:,:,1:nlev) = diff_t_entr(:,:,1:nlev)
      used = send_data ( id_diff_t_entr, diag3, Time_next, is, js, 1, mask=lmask )
      endif
      if ( id_diff_m_entr > 0 ) then
       diag3(:,:,1:nlev) = diff_m_entr(:,:,1:nlev)
      used = send_data ( id_diff_m_entr, diag3, Time_next, is, js, 1, mask=lmask )
      endif
 endif

!--- geopotential height relative to the surface on full and half levels ----

   if ( id_z_half > 0 ) then
      !--- set up local mask for fields with surface data ---
      if ( present(mask) ) then
         lmask(:,:,1)        = .true.
         lmask(:,:,2:nlev+1) = mask(:,:,1:nlev) > 0.5
      else
         lmask = .true.
      endif
      used = send_data ( id_z_half, z_half, Time_next, is, js, 1, mask=lmask )
   endif
   
   if ( id_z_full > 0 ) then
      used = send_data ( id_z_full, z_full, Time_next, is, js, 1, rmask=mask)
   endif
   
!--- zonal and meridional wind on mass grid -------

   if ( id_uwnd > 0 ) then
      used = send_data ( id_uwnd, uu, Time_next, is, js, 1, rmask=mask)
   endif
  
   if ( id_vwnd > 0 ) then
      used = send_data ( id_vwnd, vv, Time_next, is, js, 1, rmask=mask)
   endif
  
! --->h1g, 2012-08-07, dump whether stable-scheme is on, clubb_on
   if( do_clubb > 0) then
     if ( id_clubb_on > 0 ) then
        used = send_data ( id_clubb_on, clubb_on, Time_next, is, js, 1, mask=lmask )
     endif
   endif
! <---h1g, 2012-08-07
   
!-----------------------------------------------------------------------

end subroutine vert_turb_driver

!#######################################################################

subroutine vert_turb_driver_init (lonb, latb, id, jd, kd, axes, Time, &
                                  Exch_ctrl, Physics_control, &
                                  doing_edt, doing_entrain, do_clubb_in)

!-----------------------------------------------------------------------
   real, dimension(:,:), intent(in) :: lonb, latb
   integer,         intent(in) :: id, jd, kd, axes(4)
   type(exchange_control_type), intent(in) :: Exch_ctrl
   type(physics_control_type), intent(in) :: Physics_control
   type(time_type), intent(in) :: Time
   logical,         intent(out) :: doing_edt, doing_entrain

!-->h1g
   integer, optional,    intent(in)    :: do_clubb_in
!<--h1g
!-----------------------------------------------------------------------
   integer, dimension(3) :: full = (/1,2,3/), half = (/1,2,4/)
   integer :: ierr, unit, io, logunit

      if (module_is_initialized)  &
          call error_mesg  &
                   ('vert_turb_driver_init in vert_turb_driver_mod',  &
                    'attempting to call initialization twice', FATAL)

!-----------------------------------------------------------------------
!--------------- read namelist ------------------

#ifdef INTERNAL_FILE_NML
   read (input_nml_file, nml=vert_turb_driver_nml, iostat=io)
   ierr = check_nml_error(io,'vert_turb_driver_nml')
#else   
      if (file_exist('input.nml')) then
         unit = open_namelist_file (file='input.nml')
         ierr=1; do while (ierr /= 0)
            read  (unit, nml=vert_turb_driver_nml, iostat=io, end=10)
            ierr = check_nml_error (io, 'vert_turb_driver_nml')
         enddo
  10     call close_file (unit)
      endif
#endif

!---------- output namelist --------------------------------------------

      logunit = stdlog()
      if ( mpp_pe() == mpp_root_pe() ) then
           call write_version_number(version, tagname)
           write (logunit,nml=vert_turb_driver_nml)
      endif

      use_tau = Physics_control%use_tau

!     --- check namelist option ---
      if ( trim(gust_scheme) /= 'constant' .and. &
           trim(gust_scheme) /= 'beljaars' ) call error_mesg &
         ('vert_turb_driver_mod', 'invalid value for namelist '//&
          'variable GUST_SCHEME', FATAL)

      if (do_molecular_diffusion .and. do_mellor_yamada)  &
         call error_mesg ( 'vert_turb_driver_mod', 'cannot activate '//&
              'molecular diffusion with mellor_yamada', FATAL)
 
      if (do_molecular_diffusion .and. do_tke_turb)  &
         call error_mesg ( 'vert_turb_driver_mod', 'cannot activate '//&
              'molecular diffusion with tke_turb', FATAL)
 
      if (do_tke_turb .and. do_mellor_yamada)  &
         call error_mesg ( 'vert_turb_driver_mod', 'cannot activate '//&
              'tke_turb with mellor_yamada', FATAL)
 
       if (do_molecular_diffusion .and. do_edt)  &
         call error_mesg ( 'vert_turb_driver_mod', 'cannot activate '//&
           'molecular diffusion with EDT', FATAL)


!----------------------------------------------------
!   get the number of prognostic tracers
!   use later to determine prognostic vs. diagnostic
!----------------------------------------------------
      call get_number_tracers (MODEL_ATMOS, num_prog=ntp)

!-----------------------------------------------------------------------
! -->h1g, 2012-07-16
    if (present(do_clubb_in)) then
         do_clubb = do_clubb_in
    else
         do_clubb = 0
    endif
    
    doing_prog_clouds = Exch_ctrl%doing_prog_clouds

    if( do_entrain .and. do_clubb>0 ) &
         call error_mesg ( 'vert_turb_driver_mod', 'cannot activate '//&
           'both do_entrain and CLUBB', FATAL)
    nwp2 = get_tracer_index ( MODEL_ATMOS, 'wp2' )
    if ( do_clubb>0 .and. nwp2 <= ntp ) then
     ! nwp2 is a diagnostic tracer
      call error_mesg ('vert_turb_driver_mod', &
                      'wp2 is a diagnostic tracer in CLUBB', FATAL)
    endif
!<--h1g, 2012-07-16

! get tracer indices for stratiform cloud variables
          nql = get_tracer_index ( MODEL_ATMOS, 'liq_wat' )
          nqi = get_tracer_index ( MODEL_ATMOS, 'ice_wat' )
          nqa = get_tracer_index ( MODEL_ATMOS, 'cld_amt' )
       if (doing_prog_clouds) then
          if (mpp_pe() == mpp_root_pe()) &
                 write (logunit,'(a,3i4)') 'Stratiform cloud tracer indices: nql,nqi,nqa =',nql,nqi,nqa
          if (min(nql,nqi,nqa) <= 0) call error_mesg ('moist_processes', &
                         'stratiform cloud tracer(s) not found', FATAL)
          if (nql == nqi .or. nqa == nqi .or. nql == nqa) call error_mesg ('moist_processes',  &
       'tracers indices cannot be the same (i.e., nql=nqi=nqa).', FATAL)
      endif

!----------------------------------------------------------------------

      if (do_mellor_yamada) call my25_turb_init (id, jd, kd)

      if (do_tke_turb) then
        ntke = get_tracer_index ( MODEL_ATMOS, 'tke' )
        ! tke must be a diagnostic tracer
        if (ntke <= ntp) call error_mesg ('vert_turb_driver_mod', &
                    'tke can not be a prognostic tracer', FATAL)
        call tke_turb_init (lonb, latb, axes, Time, id, jd, kd)
      end if

      if (do_shallow_conv)  call shallow_conv_init (kd)

      if (do_stable_bl)     call stable_bl_turb_init ( axes, Time )

      if (do_edt)           call edt_init (lonb, latb, axes,Time,id,jd,kd)

      if (do_entrain)       call entrain_init (lonb, latb, axes,Time,id,jd,kd)
      
!-----------------------------------------------------------------------
!----- initialize diagnostic fields -----

   id_uwnd = register_diag_field ( mod_name, 'uwnd', axes(full), Time, &
        'zonal wind on mass grid', 'meters/second' ,                   &
         missing_value=missing_value    )

   id_vwnd = register_diag_field ( mod_name, 'vwnd', axes(full), Time, &
        'meridional wind on mass grid', 'meters/second' ,              &
        missing_value=missing_value    )

   id_z_full = &
   register_diag_field ( mod_name, 'z_full', axes(full), Time,    &
        'geopotential height relative to surface at full levels', &
         'meters' , missing_value=missing_value    )

   id_z_half = &
   register_diag_field ( mod_name, 'z_half', axes(half), Time,    &
        'geopotential height relative to surface at half levels', &
        'meters' , missing_value=missing_value    )

if (do_mellor_yamada .or. do_tke_turb) then

   id_tke = &
   register_diag_field ( mod_name, 'tke', axes(half), Time,      &
                        'turbulent kinetic energy',  'm2/s2'   , &
                        missing_value=missing_value               )

   id_lscale = &
   register_diag_field ( mod_name, 'lscale', axes(half), Time,    &
                        'turbulent length scale',  'm'   ,        &
                        missing_value=missing_value               )

   id_lscale_0 = &
   register_diag_field ( mod_name, 'lscale_0', axes(1:2), Time,   &
                        'master length scale',  'm'               )
endif

 if (do_edt) then
 
   id_tke = &
   register_diag_field ( mod_name, 'tke', axes(half), Time,      &
                         'turbulent kinetic energy',  'm2/s2'   , &
                         missing_value=missing_value               )
 
  end if

   id_z_pbl = &
   register_diag_field ( mod_name, 'z_pbl', axes(1:2), Time,       &
                        'depth of planetary boundary layer',  'm'  )
!--- cmip variable name ---
   id_bldep = register_cmip_diag_field_2d ( mod_name, 'bldep', Time, &
                                        'Boundary Layer Depth', 'm', &
                        standard_name = 'atmosphere_boundary_layer_thickness' )

!-->cjg: addition for new PBL depth diagnostic
   id_tref = &
   register_diag_field ( mod_name, 'tref', axes(1:2), Time,       &
                        'surface air temperature',  'K'  )
   id_qref = &
   register_diag_field ( mod_name, 'qref', axes(1:2), Time,       &
                        'surface air specific humidity',  'kg/kg'  )
   id_z_Ri_025 = &
   register_diag_field ( mod_name, 'z_Ri_025', axes(1:2), Time,       &
                        'Critical bulk Richardson height',  'm'  )
!<--cjg

!-->h1g: add RH (%) at new PBL depth,  2015-04-02 
   id_rh_Ri_025 = &
   register_diag_field ( mod_name, 'rh_Ri_025', axes(1:2), Time,       &
                        'Relative humidity at the critical bulk Richardson height',  '%'  )
!<--h1g,  2015-04-02


   id_gust = &
   register_diag_field ( mod_name, 'gust', axes(1:2), Time,        &
                        'wind gustiness in surface layer',  'm/s'  )

   id_diff_t = &
   register_diag_field ( mod_name, 'diff_t', axes(half), Time,    &
                        'vert diff coeff for temp',  'm2/s'   ,   &
                        missing_value=missing_value               )

   id_diff_m = &
   register_diag_field ( mod_name, 'diff_m', axes(half), Time,      &
                        'vert diff coeff for momentum',  'm2/s'   , &
                        missing_value=missing_value               )

! cmip6 diagnostics

   ID_evu = register_cmip_diag_field_3d ( mod_name, 'evu', Time, &
                      'Eddy Viscosity Coefficient for Momentum Variables', 'm2 s-1', &
                       standard_name='atmosphere_momentum_diffusivity', &
                       interp_method = 'conserve_order1', axis='half' )

   ID_edt = register_cmip_diag_field_3d ( mod_name, 'edt', Time, &
                      'Eddy Diffusivity Coefficient for Temperature Variable', 'm2 s-1', &
                       standard_name='atmosphere_heat_diffusivity', &
                       interp_method = 'conserve_order1', axis='half' )
                     

if (do_shallow_conv) then

   id_diff_sc = &
   register_diag_field ( mod_name, 'diff_sc', axes(half), Time,      &
                        'vert diff coeff for shallow conv', 'm2/s' , &
                        missing_value=missing_value               )
endif

if (do_stable_bl) then
  id_diff_t_stab = &
    register_diag_field ( mod_name, 'diff_t_stab', axes(half), Time,       &
                       'vert diff coeff for temp',  'm2/s',                &
                        missing_value=missing_value               )

  id_diff_m_stab = &
    register_diag_field ( mod_name, 'diff_m_stab', axes(half), Time,       &
                       'vert diff coeff for momentum',  'm2/s',            &
                       missing_value=missing_value               )

  id_stable_on = &
    register_diag_field ( mod_name, 'stable_on', axes(half), Time,       &
   'frequency of stable is on (diff_m_stab > diff_m  or diff_t_stab > diff_t)',  ' ', &
                        missing_value=missing_value               )

  id_stable_effective = &
    register_diag_field ( mod_name, 'stable_effective', axes(half), Time,       &
   'frequency of stable is effective (diff_m_stab > diff_m or diff_t_stab > diff_t and diff_m_stab > diff_min or diff_t_stab > diff_min)',  ' ', &
                        missing_value=missing_value               )

 endif


if (do_entrain) then
  id_diff_m_entr = &
    register_diag_field ( mod_name, 'diff_m_entr', axes(half), Time,        &
            'momentum vert diff coeff from entrainment module',  'm2/s',    &
                        missing_value=missing_value               )

  id_diff_t_entr = &
    register_diag_field ( mod_name, 'diff_t_entr', axes(half), Time,        &
            'heat vert diff coeff from entrainment module',  'm2/s',        &
                        missing_value=missing_value               )

 endif

! --->h1g, 2012-08-07, register id_clubb_on
if( do_clubb > 0 ) then
   id_clubb_on = &
    register_diag_field ( mod_name, 'clubb_on', axes(half), Time,       &
                       'frequency of clubb is on',  ' ',                &
                        missing_value=missing_value               )
endif
! <---h1g, 2012-08-07


!-----------------------------------------------------------------------

   doing_edt = do_edt
   doing_entrain = do_entrain
   module_is_initialized =.true.

!-----------------------------------------------------------------------

end subroutine vert_turb_driver_init


!#######################################################################

subroutine vert_turb_driver_end

!-----------------------------------------------------------------------
      if (do_mellor_yamada) call my25_turb_end
      if (do_tke_turb)      call tke_turb_end
      if (do_edt) call edt_end
      if (do_entrain) call entrain_end
      module_is_initialized =.false.

!-----------------------------------------------------------------------

end subroutine vert_turb_driver_end

!#######################################################################
! <SUBROUTINE NAME="vert_turb_driver_restart">
!
! <DESCRIPTION>
! write out restart file.
! Arguments: 
!   timestamp (optional, intent(in)) : A character string that represents the model time, 
!                                      used for writing restart. timestamp will append to
!                                      the any restart file name as a prefix. 
! </DESCRIPTION>
!
subroutine vert_turb_driver_restart(timestamp)
  character(len=*), intent(in), optional :: timestamp

   if (do_mellor_yamada) call my25_turb_restart(timestamp)
end subroutine vert_turb_driver_restart
! </SUBROUTINE> NAME="vert_turb_driver_restart"


!-->cjg: addition for new PBL depth diagnostic

subroutine bulk_Ri_height_a(t,qv,u,v,pfull,zfull,phalf,zhalf,z_Ri_025)

! This subroutine computes the height where the bulk Richardson number
! reaches 0.25. This height can be regarded an approximation to the
! PBL height. The algorithm was constructed to enable straight-forward 
! comparisons with radiosondes data.
!
! Reference: Seidel et al. 2012 (doi:10.1029/2012JD018143)
 
!-----------------------------------------------------------------------
!
!      variables
!
!      -----
!      input
!      -----
!
!      t         temperature (K)
!      qv        water vapor specific humidity (kg vapor/kg air)
!      u         zonal wind (m/s)
!      v         meridional wind (m/s)
!      zfull     height of full levels (m)
!      zhalf     height of half levels (m)
!      pfull     pressure at full levels (Pa)
!      phalf     pressure at half levels (Pa)
!
!      ------------
!      output
!      ------------
!
!      z_Ri_025  bluk Richardson height

! Input/output variables

real,            intent(in),    dimension(:,:,:) :: t,qv,u,v
real,            intent(in),    dimension(:,:,:) :: pfull, zfull
real,            intent(in),    dimension(:,:,:) :: phalf, zhalf
real,            intent(out),   dimension(:,:)   :: z_Ri_025

! Local variables

integer                                          :: i, j, k
integer                                          :: nlev, nlat, nlon
real, dimension(size(t,1),size(t,2))             :: zsurf
real, dimension(size(t,1),size(t,2),size(t,3))   :: zfull_ag
real, dimension(size(t,1),size(t,2),size(t,3)+1) :: zhalf_ag

real thetavs, thetavh, zs, us, vs, vv
real, dimension(size(t,3)) :: Ri

! Constants
real, parameter :: eps = 1.0e-8
real, parameter :: Ri_crit = 0.25
real, parameter :: p00 = 1000.0e2

!-----------------------------------------------------------------------

! Initialization
  nlev = size(t,3)
  nlat = size(t,2)
  nlon = size(t,1)
  z_Ri_025(:,:) = missing_value
       
! Compute height above surface
  zsurf(:,:) = zhalf(:,:,nlev+1)
  do k = 1, nlev
     zfull_ag(:,:,k) = zfull(:,:,k) - zsurf(:,:)
     zhalf_ag(:,:,k) = zhalf(:,:,k) - zsurf(:,:)
  end do
  zhalf_ag(:,:,nlev+1) = zhalf(:,:,nlev+1) - zsurf(:,:)
       
! Horizontal loop
  do j=1,nlat
    do i=1,nlon

!      Lowest model level properties
       thetavs = t(i,j,nlev) * (p00/pfull(i,j,nlev))**kappa  &
                             * ( 1.0 + 0.61*qv(i,j,nlev)/(1.0-qv(i,j,nlev)) )
       zs = zfull_ag(i,j,nlev)
       us = 0.0
       vs = 0.0
       Ri(nlev) = 0.0

!      Vertical upward loop
       do k = nlev-1,1,-1

         thetavh = t(i,j,k) * (p00/pfull(i,j,k))**kappa   &
                            * ( 1.0 + 0.61*qv(i,j,k)/(1.0-qv(i,j,k)) )
         vv = max( (u(i,j,k)-us)**2 + (v(i,j,k)-vs)**2, eps )
         Ri(k) = grav * (thetavh-thetavs) * (zfull_ag(i,j,k)-zs) / (thetavs*vv)
         if (Ri(k) >= Ri_crit) then
           z_Ri_025(i,j) = zfull_ag(i,j,k+1)   &
                           + (Ri_crit-Ri(k+1))/(Ri(k)-Ri(k+1))*(zfull_ag(i,j,k)-zfull_ag(i,j,k+1))
           exit
         end if

       end do

    end do
  end do

end subroutine bulk_Ri_height_a

subroutine bulk_Ri_height_b(t,qv,u,v,t_ca,q_ca,pfull,zfull,phalf,zhalf,z_Ri_025)

! This subroutine computes the height where the bulk Richardson number
! reaches 0.25. This height can be regarded an approximation to the
! PBL height. The algorithm was constructed to enable straight-forward 
! comparisons with radiosondes data.
!
! Reference: Seidel et al. 2012 (doi:10.1029/2012JD018143)
 
!-----------------------------------------------------------------------
!
!      variables
!
!      -----
!      input
!      -----
!
!      t         temperature (K)
!      qv        water vapor specific humidity (kg vapor/kg air)
!      u         zonal wind (m/s)
!      v         meridional wind (m/s)
!      t_ca      canopy air temperature (K)
!      q_ca      canopy air specific humidity (kg vapor/kg air)
!      zfull     height of full levels (m)
!      zhalf     height of half levels (m)
!      pfull     pressure at full levels (Pa)
!      phalf     pressure at half levels (Pa)
!
!      ------------
!      output
!      ------------
!
!      z_Ri_025  bluk Richardson height

! Input/output variables

real,            intent(in),    dimension(:,:,:) :: t, qv, u, v
real,            intent(in),    dimension(:,:)   :: t_ca, q_ca
real,            intent(in),    dimension(:,:,:) :: pfull, zfull
real,            intent(in),    dimension(:,:,:) :: phalf, zhalf
real,            intent(out),   dimension(:,:)   :: z_Ri_025

! Local variables

integer                                          :: i, j, k
integer                                          :: nlev, nlat, nlon
real, dimension(size(t,1),size(t,2))             :: zsurf
real, dimension(size(t,1),size(t,2),size(t,3))   :: zfull_ag
real, dimension(size(t,1),size(t,2),size(t,3)+1) :: zhalf_ag

real thetavs, thetavh, zs, us, vs, vv
real, dimension(size(t,3)+1) :: Ri, z

! Constants
real, parameter :: eps = 1.0e-8
real, parameter :: Ri_crit = 0.25
real, parameter :: p00 = 1000.0e2

!-----------------------------------------------------------------------

! Initialization
  nlev = size(t,3)
  nlat = size(t,2)
  nlon = size(t,1)
  z_Ri_025(:,:) = missing_value
       
! Compute height above surface
  zsurf(:,:) = zhalf(:,:,nlev+1)
  do k = 1, nlev
     zfull_ag(:,:,k) = zfull(:,:,k) - zsurf(:,:)
     zhalf_ag(:,:,k) = zhalf(:,:,k) - zsurf(:,:)
  end do
  zhalf_ag(:,:,nlev+1) = zhalf(:,:,nlev+1) - zsurf(:,:)
       
! Horizontal loop
  do j=1,nlat
    do i=1,nlon

!      Extended height
       z(1:nlev) = zfull_ag(i,j,:)
       z(nlev+1) = 2.0         ! Assumed reference (canopy) air height

!      Canopy level properties
       thetavs = t_ca(i,j) * (p00/phalf(i,j,nlev+1))**kappa  &
                           * ( 1.0 + 0.61*q_ca(i,j)/(1.0-q_ca(i,j)) )
       zs = z(nlev+1)
       us = 0.0
       vs = 0.0
       Ri(nlev+1) = 0.0

!      Vertical upward loop
       do k = nlev,1,-1

         thetavh = t(i,j,k) * (p00/pfull(i,j,k))**kappa   &
                            * ( 1.0 + 0.61*qv(i,j,k)/(1.0-qv(i,j,k)) )
         vv = max( (u(i,j,k)-us)**2 + (v(i,j,k)-vs)**2, eps )
         Ri(k) = grav * (thetavh-thetavs) * (z(k)-zs) / (thetavs*vv)
         if (Ri(k) >= Ri_crit) then
           z_Ri_025(i,j) = z(k+1) + (Ri_crit-Ri(k+1))/(Ri(k)-Ri(k+1))*(z(k)-z(k+1))
           exit
         end if

       end do

    end do
  end do

end subroutine bulk_Ri_height_b
!<--cjg

end module vert_turb_driver_mod

